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We present a Bayesian dynamical inference method for characterizing cardiorespiratory (CR) 
dynamics in humans by inverse modelling from blood pressure time-series data. This new method is 
applicable to a broad range of stochastic dynamical models, and can be implemented without severe 
computational demands. A simple nonlinear dynamical model is found that describes a measured 
blood pressure time-series in the primary frequency band of the CR dynamics. The accuracy of 
the method is investigated using surrogate data with parameters close to the parameters inferred 
in the experiment. The connection of the inferred model to a well-known beat-to-beat model of the 
baroreflex is discussed. 
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I. INTRODUCTION 

Model identification is an important method used in 
both fundamental and applied research 0, HIE El IE IE 
013 on the human cardiovascular system (CVS). Because 
of the complexity of CVS dynamics and the multiplicity 
of its mechanism, it is inherently difficult or impossible 
to isolate and study individual response mechanisms in 
the intact organism [9j . In such cases mathematical mod- 
els of cardiovascular control that are consistent with the 
experimental data can provide valuable insights 0, ^| . 
Altered dynamics of the cardiovascular system is asso- 
ciated with a range of cardiovascular diseases and with 
increased mortality, and it is hoped that dynamical met- 
rics will provide new means of evaluating autonomic ac- 
tivity, and eventually form the basis for diagnostic tests 
for many conditions (see e.g. 0,0,0]). 

Despite the fact that most cardiovascular controls are 
demonstrably nonlinear |H0 0J30 l2l| and 
are perturbed by stochastic inputs |c |j2a.j 2a | . oversim- 
plified assumptions of model linearity fl1.l3ll3.IMlri lililT(| 

and / or determinism jToL Ht| are often made in order to 
make some progress in cardiovascular system identifica- 
tion. Such choices are often influenced more by the avail- 
ability of certain statistical tools and methodologies than 
by biophysical or medical considerations. It is very de- 
sirable to develop reliable methods of system identifica- 
tion that do not have such limitations, and are capa- 
ble of treating more realistic models. Such models could 
be used to relate difficult-to-access parameters to non- 
invasively-measured data |Tl|. 

While a number of numerical schemes have been pro- 
posed recently to deal with different aspects of the in- 
verse problem using linear approximations Q, |E 13 HE] , 
or semi-quantitative estimations of either the strength 
of some of the nonlinear terms [2(j or the directional- 
ity of coupling |5E l2c| . inverse cardiovascular inverse 
problems remain difficult because of the complexity and 



nonlincarity of the cardiovascular interactions, as well 
as the stochasticity of many dynamical inputs to the 
system. The problem of nonlinear cardiovascular sys- 
tem identification has been addressed in a number of 
publications 0, IE E3- Nonlinearities generally require 
the use of more complex and involved numerical tech- 
niques HE EE El HI HI EE ED> while the presence 
of dynamical noise in continuous systems can introduced 
systematic errors in the estimation of the model param- 
eters EE 13- Analogous difficulties arise in a broad 
range of problems in many scientific disciplines, including 
problems in lasers EH and molecular motors EE] , m epi- 
demiology EE] , an d m coupled matter-radiation systems 
in astrophysics Ell- An obstacle to progress in these 
fields is the lack of general methods of dynamical infer- 
ence for stochastic nonlinear systems. Accordingly the 
methods described in this paper should be of broad in- 
terdisciplinary interest. 

In this paper we introduce a novel method for the 
analysis of cardiorespiratory dynamics within a nonlin- 
ear Bayesian framework for the inference of stochastic 
dynamical systems 34] ■ This method is applied to the 
analysis of a univariate blood pressure time-series where 
a simple nonlinear dyna mical model based on coupled 
nonlinear oscillators 0, 142. |43| . is found that describes 
time-series data in the relevant frequency range. The ac- 
curacy of the method is investigated using surrogate data 
with parameters close to the parameters inferred in the 
experiment, and the connection of this model to a well- 
known beat-to-beat model of the baroreflex is discussed. 



II. THE METHODOLGY 

In the methodological framework presented here there 
are three essential steps: (i) the input data is prepared, 
(ii) a parameterized class of models chosen that describes 
the data, and (iii) the parameters of this model are in- 
ferred from the time-series data. 
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A. Data 

Here we worked with a particular recording of the cen- 
tral venous blood pressure, a sample of which is shown 
in the Fig. [TJa). A feature this blood pressure (BP) 
time-series is the presence of the two oscillatory com- 
ponents at frequencies approximately f r rs 0.2 Hz and 
f c ~ 1.7 Hz corresponding to the respiratory and cardiac 
oscillations. It can also be clearly seen from the spec- 
tra that the nonlinear terms including terms of nonlinear 
cardiorespiratory interaction (corresponding to the side 
peaks) are very visible in this sample. We note that the 
relative intensity and position of the cardiac and respi- 
ratory components vary strongly from sample to sample 
with average frequency of the respiration being around 
0.3 Hz and of the heart beat being around 1.1 Hz. 

In preparing cardiovascular data for model identifica- 
tion one has to bear in mind that the CVS power spectra 
reflect a large variety of complex cardiovascular interac- 
tions seen as peaks and other features in a very broad 
frequency range [22|, |6(J, uH], |63] • In order to make sense 
of these multi-scale phenomena parametric modelling is 
usually restricted to a specific part of the power spec- 
trum. It is clear that in modelling the cardiorespiratory 
interaction the frequency range of modelling must include 
at least main harmonics of cardiac and respiratory oscil- 
lations f c and f r and their combinational frequencies. 
Moreover, as was pointed out already by Womersley (see 
e.g. [6l| cf. also with |g4|) locally measured blood pres- 
sure signals resembles a steady-state oscillations and the 
sum of the first three harmonics contains more then 70% 
of the total signal variance. Therefore, it is desirable that 
at least three harmonics (see also discussion below) of the 
basic frequencies of the respiratory and cardiac oscilla- 
tions are included into the frequency range of modelling. 



B. Models 

When one considers modelling the cardiovascular sys- 
tem, one usually envisions constructing a model based 
on biophysical principles that is capable of generating 
solutions that reproduce, to some degree, the data: the 
forward modelling problem (see e.g. 0, 0, 
One may also consider the inverse modelling problem, in 
which models are built that describe measured data (see 
e -g- fltl3ll4l.l5l.l32j). Both approaches have proven useful 
in the context of the cardiovascular research with forward 
approach providing valuable insight into the system and 
its causal relationships, and the inverse approach provid- 
ing a useful means of intelligent patient monitoring of 
cardiovascular function. 

As a third alternative one may try to bridge the two ap- 
proaches by building a model that accurately reproduces 
the experimental observations while at the same time is 
based on the biophysical principles of circulation. In such 
a case the form of the mathematical model is taken from 
biophysical principles, with its component parts corre- 



sponding to a greater or lesser degree to specific biophys- 
ical mechanisms, while the values of some or all of the pa- 
rameters of the mathematical model are inferred directly 
from from the data. In such a case it is to be hoped that 
information with direct biophysical significance, and not 
mere mathematical or statistical characterizations can be 
inferred from the data. 

Many studies have been carried out to explore the 
physiological mechanisms underlying cardiorespiratory 
interactions HI. 1501 l5l| . The most important ones are 
the modulation of cardiac filling pressure by respira- 
tory movements [52| , the direct respiratory modulation of 
parasympathetic and sympathetic neural activity in the 
brain stem |53j , and the respiratory modulation of the 
baroreceptor feedback control [54J. A common feature 
that these mechanisms is that they are nonlinear, have 
a dynamical (or memory) com ponent, and are subject to 
exogenous fluctuations 0, liHl42l I5M 15(1 l57l l58| . 

A simple beat-to-beat model describing the cardio- 
respiratory systems DeBoer et al. 0, |5!| . The DeBoer 
model has further been elaborated recently in [ToL ITTL f23| . 
Insight into cardio-respiratory dynamics can also be 
gained through inverse modelling where the cardiac and 
respiratory cycles are modelled in terms of coupled non- 
linear oscillators 0, |2^, E2> CU EH • In this approach 
spectral and synchronization features observed in the 
time-series data are interpreted physiologically, and re- 
lated to the model parameters [Tjj- However, the iden- 
tification of the model parameters could not be inferred 
directly from the time-series data. Instead an extensive 
computer simulations have been employed to establish 
realistic values for the model parameters p2|. 

The simplest model that could reproduce steady-state 
oscillations of the blood pressure signal at two fundamen- 
tal frequencies is a system of two coupled limit cycles 
on a plane. According to the results of the Poincare- 
Bcndixson theory of planar dynamical systems for a sys- 
tem to have limit cycle in a simply connected region the 
divergence of the vector field must change sign in this 
region (see e.g. [HJ). Therefore, we conclude, that the 
simplest system that can reproduced the discussed fea- 
tures of the BP signal is planar systems with limit cycles 
which vector field contains polynomials of the order 3. 
Accordingly, we model the time-series data as a system of 
two coupled oscillators with vector fields including non- 
linearities (including nonlinearities in coupling terms) up 
to the 3-rd order in the form 

N 2 

x r = a x x r + hy r , y r = y^a^x, y) +y"Vij$j,(l) 
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x c = a 2 x c + b 2 yc, Vc = ^ 0i<f>i(x, y) + ^°'2jtj, (2) 

i=i j=i 

&(*)> = 0, te(*)0(*')> = M(*-*')- 

Here noise matrix a mixes zero-mean white Gaussian 
noises £j(t), which is related to the diffusion matrix 
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FIG. 1: Data from record 24 time series of the MGH/MF Waveform Database available at www.physionet.org. (a) Original 
time series of the central venous blood pressure and (b) Power spectrum of original time series, (c) Respiratory component 
produced by filtering the blood pressure time series with a 0.06 Hz, Order 2, zero-phase, high-pass Butterworth filter and 
a 0.6 Hz, order f2, zero phase low-pass Butterworth filter, and (d) the power spectrum of the respiratory component, (e) 
Cardiac component produced by filtering the blood pressure time series with a 0.8-3.0 Hz Hz, order 8, zero-phase, band-pass 
Butterworth filter and (d) the power spectrum of the cardiac component. The chosen frequency range of of the components 
selecting according to the discussion in the text. 



D = ao T . And base functions are chosen in the form 

jfi 2222 3 

L > T - ' c? Vr i Vc i 3 % c ' Vr ' Vc ' ^rUr 3 %cVc 1 %r ' 
322 2233 2 2l fn\ 

The restrictions imposed on the rhs of the equations for 
x r and x c in and J21 are determined mainly by the 
fact that we have to infer four hidden dynamical vari- 
ables using univariate time-series data (see next section 
for further details). The parametric representation 
and J2J covers a wide range of models with limit cycles 
in the plane. In particular, with a special choice of the 
model parameters it describes van der Pol or FitzHugh- 
Nagumo oscillator systems that are very popular in the 
context of the cardiovascular modelling. Furthermore, 
the choice of the parametric model in the form and 
(0) allows one to relate them to physiological parameters 



characterizing autonomous nervous system (see section 
I VII for the discussion). See 19] for the alternative choice 
and corresponding physiological reasoning. 

C. Parameters 

Following the logic of the inverse modelling ap- 
proach, we must then identify the parameters M. — 
{a, b, a, /3, D} of the model © that reproduce the 
dynamical and spectral features of the BP signal shown 
in the Fig. 1. Terms representing nonlinear cardio- 
respiratory interactions are described by the last three 
base functions in The correspondence of these terms 
to the experimentally observed combinational frequen- 
cies in the BP signal is summarized in the Fig. [21 It 
can be seen from the figure that the same combinational 
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FIG. 2: Summary of the main harmonics of the cardiac and 
respiratory components observed in the BP signal. The cor- 
respondence between the nonlinear terms of the model Q, 
J2J) and the frequencies observed in the time-series data are 
shown by arrows. 



frequencies correspond to the nonlinear coupling terms 
in both limit cycle systems in the model therefore a non- 
linear time-series analysis is a requirement for the iden- 
tification of such a model. 

Here we show how the task of identifying a model of the 
cardiovascular system based on coupled nonlinear oscil- 
lators can be performed systematically within a Bayesian 
statistical framework. In this approach parameters of the 
model can be inferred directly from the time-series data. 
We also discuss briefly how the links between the coupled 
oscillator model and beat-to-beat models can in principle 
be established. 



III. BAYESIAN INFERENCE OF STOCHASTIC 
NONLINEAR DYNAMICAL MODELS 

Details of our new Bayesian technique can be found 
elsewhere |34j . Here we give a brief description of the 
main steps of the algorithm. 

Stochastic nonlinear dynamical models of the type 
PJl can be expressed as a multi-dimensional nonlinear 
Langevin equation 



i(i) = f(x)+e(t)=f(x)+o£(t), 



(4) 



where e(t) is an additive stationary white, Gaussian vec- 
tor noise process characterized by 



(£(*)> = 0, mZ T (t'))=V5(t-t'), 



(5) 



where D is a diffusion matrix. 

It is assumed that the trajectory x(t) of this system is 
observed at sequential time instants {if.; k = 0, 1, . . . , K} 
and a series S = {s k = s(t k )} thus obtained is related to 
the (unknown) "true" system states X = {x k = x{tk)} 
through some conditional PDF p Q (S\X). 

An a priori expert knowledge about the model param- 
eters is summarized in so-called prior PDF p pl (M). In 



our case we chose prior in the form of zero-mean Gaussian 
distribution for model parameters and uniform distribu- 
tions for the coefficients of diffusion matrix. 

If experimental time-series data S are available they 
can be used to improve the estimation of the model pa- 
rameters. The improved knowledge of the models pa- 
rameters is summarized in the posterior conditional PDF 
Ppost(M\S) , which is related to prior via Bayes' theorem: 



p post (M\S) 



£(S\M) Ppr (M) 
J e(S\M) Ppi (M)dM' 



(G) 



Here £(S\A4), usually termed the likelihood, is the con- 
ditional PDF that relates measurements S to the dy- 
namical model. The determinant on the right hand 
side of the equation 10 is merely a normalization fac- 
tor. In practice, © can be applied iteratively using a 
sequence of data blocks S,S', etc. The posterior com- 
puted from block S serves as the prior for the next block 
S' , etc. For a sufficiently large number of observations, 
Ppast(-M\S, S' , . . .) is sharply peaked at a certain most 
probable model M. = M* . 

The main efforts in the research on stochastic nonlinear 
dynamical inference are focused on constructing the like- 
lihood function, which compensates noise induced errors, 
and on introducing efficient algorithms of optimization of 
the likelihood function and integration of the normaliza- 
tion factor (cf. [33. l35l Iririlp . 

In our earlier work |34j a novel technique of nonlinear 
dynamical inference of stochastic systems was introduced 
that solves both problems. To avoid extensive numerical 
methods of optimization of the likelihood function and 
integration of the normalization factor we suggested to 
parameterize the vector field of (@J in the form 



f(x)=U(x)c = f(x;c), 



(7) 



where U(x) is an TV x M matrix of suitably chosen basis 
functions {U nm (x); n = 1 : iV, m = 1 : M}, and c is an 
A/-dimensional coefficient vector. An important feature 
of J7J) is that, while possibly highly nonlinear in x, f (x; c) 
is strictly linear in c. 

The computation of the likelihood function can be cast 
in the form of a path integral over the random trajecto- 
ries of the system [67], • Using the uniform sampling 
scheme introduced above we can write the logarithm of 
the likelihood function in the following form for suffi- 
ciently small time step h (cf. [67ll68|'): 



K-l 



2 logl(y|M) - IndetD + A ]T [v(y*)c (8) 



K 



k=0 



+ (y k - U k c) D (y k - U fe c))j + N ln(2nh), 

which relates the dynamical variables x(t) of the system 
l fi| to the observations s(t). Here we introduce the fol- 
lowing notations U fc = U(y fe ), y fe = h^ 1 (y k+1 - y k ) and 
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vector v(x) with components 



N 



•w = E 



dU nm {x) 



m = 1 : M. 



The vector elements {c m } and the matrix elements 
{-Dnn'} together constitute a set = {c, D} of unknown 
parameters to be inferred from the measurements S. 

Choosing the prior PDF in the form of Gaussian dis- 
tribution 



Ppr(-M) = 



'det(E pr 1 ) 



1 



(2tt) m 6XP V 2^° ~ c P r ) TS P rl ( c ~ c P r ) 



(9) 

and substituting p pl (A4) and the likelihood t(S\M) 
into © yields the posterior p p0 st{M\S) — const x 
exp[—L(M\S)}, where 

L(M\S) = L s (c,D) = i Ps (D)-c T w s (D) + ic T H 5 (D)c 

(10) 

Here, use was made of the definitions 



K-l 



p s {t>) = hJ2^l D" 1 s fe + K ln(det D), 

k=0 

K-l 

w s (D) = S- 1 c pr + / l ^ 



fc=0 



v(s fe ) 



K-l 



S s (D)=S- 1 +/ l ^U^D- 1 ^ 



(11) 
(12) 
(13) 



fe=0 



The mean values of c and D in the posterior distribu- 
tion give the best estimates for the model parameters for 
a given block of data S of length K and provide a global 
minimum to L s (c, D). We handle this optimization prob- 
lem in the following way. Assume for the moment that c 
is known in (|10fl . Then the posterior distribution over D 
has a mean D post = ©s(c) that provides a minimum to 

S s (c, D) with respect to D = D T . Its matrix elements 
are 



IV. ESTIMATION OF PARAMETERS OF 
CARDIORESPIRATORY INTERACTION FROM 
UNIVARIATE TIME-SERIES DATA (REAL 
DATA) 

In order to apply algorithm l|12|) - (|15|) for the identifi- 
cation of the model of nonlinear cardio-respiratory dy- 
namics , (J2J from the univariate BP time-series of the 
type shown in Fig. 1(a) we have to extract time-series 
data corresponding to the four dynamical variable in the 
model. Accordingly we divide the total spectrum into 
a low- frequency respiratory component s r (t) and high- 
frequency cardiac component s c (t) as is shown in Fig. 
1(c) and (e). 

A discussion of the physiological relevance of this spec- 
tral separation can be found in d E3). However, it 
is perfectly correct to consider this separation a mathe- 
matical ansatz. Physiological considerations need only 
come into play at the point where we attempt place 
a specific biophysical interpretation of the model ele- 
ments. The parameters of the filters (see Fig. 1 caption) 
were chosento preserve the 2-nd and 3-rd harmonics of 
these signals. Then two x r (t) and x c (t) dynamical vari- 
ables of the model QJ, J2J) can be identified with intro- 
duced above two-dimensional time-series of observations 
s(t) = {s r (t), s c (t)}. The remaining two dynamical vari- 
ables y(t) = {y r (t),y c (t)} can be related to the observa- 
tions {s(tfc)} as follows 



b n y„(t k ) 



s n {tk + h) - s n (t k ) 



+ a n s n (t k ), (16) 



where n = r,c. The relation 116(1 is a special form of em- 
bedding that allows one to infer a wider class of dynami- 
cal models of the cardiorespiratory interactions including 
models in the form of FitzHugh-Nagumo oscillators. It is 
clear now that we have introduced the restrictions on the 
form of the rhs of the first equations in , J2J to reduce 
the number of parameters of embedding that have to be 
selected to minimize the cost (fTX>|) and provide the best 
fit to the measured time series {s(tfe)}. The correspond- 
ing simplified model of the nonlinear interaction between 
the cardiac and respiratory limit cycles can now be writ- 
ten in the form corresponding to parametrization (JJJ as 
follows 
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K-l 



©r (c) = jfl^ [h - U(s fe ) cj n [s k - U(y fe ) c 

(14) 

Alternatively, assume next that D is known, and note 
from p0(l that in this case the posterior distribution over 
c is Gaussian. Its covariance is given by 3 S (D) and the 
mean Cp OSt minimizes L s (c, D) with respect to c 



"post 



H 5 - 1 (D)w s (D). 



(15) 



We repeat this two-step optimization procedure itera- 
tively, starting from some prior values c pr and £ pr . 



y = U(s,y)c + £(i), 



(17) 



where £(t) is a two-dimensional Gaussian white noise 
with correlation matrix D, and the matrix U will have 
the following block structure 



U = 





'0i o 




'02 




'0s 






0! _ 




2 _ 




B _ 





(18) 



Here B — 22 diagonal blocks of size 2x2 formed by the 
basis functions given in and the vector of unknown 
parameters c has the length M = 2B. 

Finally, the model (fTTjl . (fT%|) has to be inferred using 
method described in the previous section. The compari- 
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FIG. 3: (a) Time series of the cardiac oscillations x c (t n ) = 
Se(tn) in arbitrary units (black line) obtained from central 
venous blood pressure. The sampling rate was 90 Hz after 
resampling of the original signal. Inferred time series of the 
cardiac oscillator (green line), (b) Power spectrum of the 
cardiac oscillations obtained from the real data (black line). 
Power spectrum of the inferred oscillations (green line), (c) 
Limit cycle of the cardiac oscillations (x c {n), y c (n) obtained 
from real data as described in the text (black line). Limit 
cycle of inferred oscillations (green line). 



FIG. 4: (a) Time series of the respiratory oscillations 
Xritn) = Srifn) in arbitrary units (black line) obtained from 
central venous blood pressure. The sampling rate was 90 Hz 
after resampling of the original signal. Inferred time series of 
the respiratory oscillator (green line), (b) Power spectrum of 
the respiratory oscillations obtained from the real data (black 
line). Power spectrum of the inferred oscillations (green line), 
(c) Limit cycle of the respiratory oscillations (x c (n),y c (n) ob- 
tained from real data as described in the text (black line). 
Limit cycle of inferred oscillations (green line). 



son between the time series of the inferred and actual car- 
diac oscillations is shown in Fig. |31 Similar results are ob- 
tained for the respiratory oscillator as shown in the Fig. 
0] In particular, the parameters of the nonlinear cou- 
pling and of the noise intensity of the cardiac oscillations 
have been estimated to have the following values 020 = 
2.2,021 = 0.27, /3 22 = -8.67, and (g(t)) = 8.13. While 
parameters of the coupling of respiratory oscillations to 
the cardiac oscillations were estimated to have the follow- 
ing values a2o = 0.12, a 2 i = 0.048, ct22 = —0.066, and 
D\\ = 0.18 . Consistent with expectations, in all experi- 
ments the parameters of the nonlinear coupling are more 
then one order of magnitude higher for the cardiac os- 
cillations as compared to their values for the respiratory 
oscillations reflecting the fact that respiration strongly 
modulates cardiac oscillations, while the opposite effect 



of the cardiac oscillations on respiration is weak. 

We have shown that it possible using these methods to 
simultaneously infer the coupling strengths and noise pa- 
rameters nonlinear cardiorespiratory dynamics directly 
from a non-invasively measured time series. We view 
this demonstration of principle as a first step towards the 
practical use this technique for cardio-respiratory mod- 
elling and and in clinical applications. A number of very 
important physiological and mathematical issues arise in 
relation to the application of this new technique to spe- 
cific problems. We hope to address many of these in 
future publications. In what follow here we specifically 
consider the problem of estimating of the accuracy of the 
method, and begin the discussion of the connection be- 
tween inferred parameters and the indexes of autonomous 
cardiovascular controls. 
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FIG. 5: (a) The velocity of the respiratory component of 
oscillations of the original signal y r (t) (green line) is com- 
pared with the signal y r {tk) (black dashed line) obtained 
as a result of filtration of s(t) followed by the embedding 
biy r (tk) = (s r (tk + h) — s r (tk))/h + a2S r (tk)- (b) Power spec- 
tra of the original velocity of the respiratory component y r (tk) 
(green line) is shown in comparison with the power spectrum 
of the recovered signal y r (tk) (black dashed line). 



V. VALIDATION OF THE METHOD USING 
SURROGATE TIME-SERIES DATA 

It is desirable to check performance of the method on 
surrogate time-series data obtained by numerically sim- 
ulating the model QJ, ©with the parameters measured 
with the CVS data. 

To this end we consider a surrogate signal x(t) = 
x r (t) + x c (t) as a time-series data input s(t) for the in- 
ference. Here x r (t), x c (t) are obtained using numerical 
simulations of the model QJ, J3J) with the parameters 
taken from the inference of the experimental BP signal 
described in the previous section. 

First we verify that the decomposition of the input 
signal s(t) into low- frequency s r and high-frequency s c 
harmonics using two band-pass Butterworth filters and 
subsequent application of the embedding procedure l|16|) 
allows one to reconstruct the original signal. In the Fig. 
we compare the velocity of the respiratory component 
of the original signal y r {t) with the the reconstructed 
velocity y r (t) . Similar results are obtained for the recon- 
struction of the high-frequency component. We notice in 
particular that the noise introduced by embedding can 
be neglected since it is more then order of magnitude 
smaller then the dynamical noise in the signal. 

Now we can apply inference procedure described in the 
previous section to estimate nonlinear coupling param- 
eters of the model from the univariate surrogate time- 
series data. The results of the estimation are summa- 
rized in the Table [Q It can be seen from the table that 
the method allows one to estimate correct order of the 
absolute values of the nonlinear coupling parameter. For 
some parameter the accuracy of the estimation is much 
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FIG. 6: (a) Surrogate time series of the respiratory oscil- 
lations Xritn) in arbitrary units (black line) obtained from 
model Q, Inferred time series of the cardiac oscillator 
(green line), (b) Power spectrum of the surrogate respiratory 
oscillations (black line) . Power spectrum of the inferred oscil- 
lations (green dashed line), (c) Limit cycle of the surrogate 
respiratory oscillations (x c (n), y c (n) (black line). Limit cycle 
of inferred oscillations (green dashed line). 



better, but in practice the correct values are not know. 
Therefore we conclude that the accuracy of the estima- 
tion of the absolute values of parameters of coupling of 
two limit cycle systems from univariate time-series data 
is within the order of magnitude. 

Similar results are obtained for the estimation of other 
parameters of the model. Using values of the model pa- 
rameters estimated from the univariate surrogate data 
one can reconstruct very closely the dynamical and spec- 
tral features of the original system as shown in the Fig. 
[H] The largest errors of estimations are obtained for 
the values of the noise intensity as shown in two last 
columns of the Tabled] This result can be easily under- 
stood taking into account that filtration of the signals has 
the strongest effect on the noise spectrum of the system. 
However, the filter-induced errors are systematic and can 
be corrected using test with surrogate data. 

The main source of error is related to the spectral de- 
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0.048 


0.27 


-0.066 


-8.67 


0.18 


8.13 


0.18 


6.32 


0.011 


0.49 


0.053 


6.03 


0.017 


3.44 


51.2% 


186.8% 


75.9% 


102.7% 


27.9% 


30.6% 


90.8% 


57.7% 



TABLE I: Absolute values of the coefficients of nonlinear 
cardiorespiratory interactions corresponding to the last three 
base functions {x r x c , x^.x c , x r x^}. Coefficients {at} cor- 
respond to the respiration coupling to cardiac rhythm. Coef- 
ficients {Pi} correspond to the cardiac oscillation coupling to 
respiration. For each set of coefficients the actual values (top 
row) are compared with the mean inferred values obtained 
from 100 blocks of time-series data x(t) — x r (t) + x c (t) with 
50000 points in each block and sampling time 0.02 sec (middle 
row). The error of estimation is shown in the bottom line. 



«20 


020 


«21 


021 


«22 


022 


Da 


£>22 


0.12 


2.20 


0.048 


0.27 


-0.066 


-8.67 


0.18 


8.13 


0.12 


2.41 


0.048 


0.28 


-0.070 


-8.61 


0.18 


8.14 


2.9% 


9.3% 


1.8% 


5.6% 


5.2% 


0.7% 


0.2% 


0.2% 



TABLE II: Absolute values of the coefficients of nonlinear 
cardiorespiratory interactions corresponding to the last three 
base functions 10 {x r x c , x^x c , x r x^.}. Coefficients {at} cor- 
respond to the respiration coupling to cardiac rhythm. Coef- 
ficients {0i} correspond to the cardiac oscillation coupling to 
respiration. For each set of coefficients the actual values (top 
row) are compared with the mean inferred values obtained 
from 100 blocks of time-series data {x r (t),x c (t),y r (t),y c (t)} 
with 160000 points in each block and sampling time 0.01 sec 
(middle row) . The error of estimation is shown in the bottom 
line. 

composition of the univariate data and it is therefore sys- 
tematic. To illustrate this point we use the original sur- 
rogate time-series data {x r (t),x e (t),y r (t),y c (i)} for two 
coupled oscillators to infer parameters of the model QJ, 
(J2J. The results of inference of the coupling parameters 
are shown in the Tab. [H] It can be seen that the values 
of the parameters can be estimated with relative error 
better then 10%. In particular, the relative error of es- 
timation of the noise intensity is now better then 4%. 
The accuracy of the estimation can be further improved 
by increasing the total time of observation of the system 
dynamics as explained in j34[. 

These results should be compared with semi- 
quantitative estimations of either relative strength of 
some of the nonlinear terms [2(| or the directionality 
of coupling [23, from bivariate time-series data. It 
becomes clear that our algorithm provides an alterna- 
tive effective approach to the solution of the problem of 
analysis of cardiovascular coupling. In particular, the re- 
sults of this section validate the application of the method 
to the experimental time-series cardiovascular data and 
demonstrate that the it is indeed possible to estimate si- 
multaneously the strength, directionality and the noise of 
nonlinear cardiorespiratory coupling form the univariate 
blood pressure signal. The accuracy of the estimation is 
within the order of the magnitude. The main source of 
errors is the decomposition of the univariate signal into 



two oscillatory components and it is, therefore, system- 
atic. Using this fact one can introduce systematic cor- 
rections to improve the results of the estimation of the 
parameters of the CR interaction from the experimen- 
tally measured BP signal. 

VI. DISCUSSION 

It is important to establish a relationship between the 
model parameters and physiological parameters of the 
cardiovascular system. A beat-to-beat model describing 
the relationships between blood pressures and respiration 
in simple, b ut phys iologically meaningful terms is the De- 
Boer model |43, |S3 • While the DeBoer model cannot de- 
scribe the dynamics within one heartbeat it does incorpo- 
rate several well-known physiological laws of the cardio- 
respiratory system. More recent extensions and modifi- 
cations of the DeBoer model have appeared [ToL TTTL |23| . 
The problem of inverse modelling was not addressed in 
this earlier work. It is therefore very desirable to con- 
nect the approach presented here with such beat-to-beat 
models. 

The DeBoer model describes the beat-to-beat evolu- 
tion of the state variables shown in the Fig. [7| (a): sys- 
tolic pressure (S), diastolic pressure (D), RR intervals 
(J), and arterial decay time (T = R x C = peripheral 
resistance x arterial compliance). Following a brief ac- 
count of the DeBoer model given in [lj| and neglecting 
for the sake of simplicity the variation of the peripheral 
resistance we can write a set of corresponding difference 
equations in the form 

Di = S i _iexp[(-2/3)/ i _i/T], (19) 
S, = A+7-fi-i + Ci+^sin(27r/t), (20) 
h = G v S[_ Tv + G P F(S', rp) + C 2 , (21) 

Here C±, C%, and C3 are constants and the sigmoidal 
nature of the baroreceptor sensitivity is accounted for by 
defining an effective Systolic pressure (S') [48| 

Si = S , o + 18arctan[(S'-5o)/18]. (22) 

The first equation (|19fl follows from the Windkessel 
model of the circulation, while the second equation l|20|) 
expresses contractile properties of the myocardium in ac- 
cordance with Starling's law that takes into account me- 
chanical effect of the circulation on the BP (see and 
e.g. [13 )• The last equation (|21|l includes explicitly two 
mechanisms of the cardiovascular control defined by their 
respective gain (G) and delay (t): (i) fast vagal control of 
the heart rate G v Sl_ r ^, (ii) slower /3-sympathetic control 
of the heart rate GpF(S' ,rp). Here F(S',t) is a linear 
weighted sum of the form 

M 

F (S\ t) = a kS[- T+k = (S'i_ T _ 2 + 2S' i _ T _ 1 ) 

k=-M 

+3Si_ T + 2S' i _ T+1 + S' t _ T+2 )/9 
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FIG. 7: (a) The BP signal in the frequency range of car- 
diac oscillations (black line). Systolic pressure (S n ), diastolic 
pressure (D n ), RR intervals {I n ), and arterial decay time 
(T = RC =) are shown for the n-th heartbeat, (b) Com- 
parison of the BP signal (thing black line) with the approx- 
imation adopted in the DeBoer model (dashed line) and the 
approximation by the FitzHugh-Nagumo model (dotted line) . 
Vertical scale has arbitrary units. 



Further we assume for simplicity that the pressure oscil- 
lations do not deviate far away from the working point 
S in IHl), i.e. S' « S. 

To establish the connection between DeBoer (|19|l - <|21[) 
model and the model Q, J2J introduced in this paper 
we note that the equations DeBoer model is a piece- wise 
approximation of the actual BP signal. In particular, it 
describes the BP signal as an exponential decay during 
2/3 part of the RR interval and linear increase during 
1/3 of the /„ as show in the Fig. [7\ (b). We note also 
that the model of the cardiac oscillations (0) resembles a 
model of the FitzHugh-Nagumo (FHN) system 



x = e(y - fix), 

y =ay + -fy 2 + Sy 3 -x + C, 



(23) 



where we have neglected for a moment the cardiorespira- 
tory interaction. The approximation of the BP signal by 
the output of the FHN system is also shown in the Fig. 




FIG. 8: Time evolution of the dynamical variables x (solid 
line) and y (dashed line) of the FHN system with the following 
parameters: e = 0.01, (3 = -0.05, C = -0.125, a = 0.5, 
7=1,5 = -1. 

0(b). It can be seen already from a comparison between 
two approximations that there is a close connections be- 
tween DeBoer model and model of coupled oscillators 
considered in this paper. This can be further illustrated 
by noticing that for small e the limit cycle in the FHN 
system consists of fast motion with practically constant 
values of y, when x jumps between negative and positive 
values, and slow motion, when x changes very little (see 
Fig.JSJ. Assuming the constant value of x at the top \a+\ 
and at the bottom — jot j of the dashed curve that cor- 
respond to the slow motion along the limit cycle we can 
integrate the first equation in (|23|l to obtain 



, (S^-la-De-^ + l. 
X ^>-S ( Dn + \a + \)e-Pt-\a A 



for < t < |J n ; 
for %I n < t < I n . 



This solution closely resembles eqs. (|19|l and (12011 of the 
DeBoer model 

It can be seen even from this simplified discussion that 
the parameters of the model 0J, J2J found in the present 
paper can be related directly to the physiological param- 
eters of the autonomous control of circulation. Further- 
more, this discussion suggests that it should be possible 
at least in principle to bridge inverse and forward mod- 
elling and to infer parameters of the autonomous ner- 
vous control of the cardiovascular system directly from 
the time-series data. 

We emphasize, however, that the obtained results is 
only the first step in this direction. In particular, the 
DeBoer model itself has to be modified in various ways, 
including more realistic functional form of the feedback 
terms and specifically taking into account the fact that 
the baroreflex control is a closed loop [TsL l2lj . In fact 
it was shown that a multi-compartment closed-loop 
model of the cardiovascular responses can simulate well 
the experimentally observed variations in the time-series. 
On the other hand, this comparison suggests that the in- 
ference scheme used in this paper has to be modified in 
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a various ways to facilitate convergence and guarantee 
deeper physiological meaning of the model parameters 
as will be discussed in more details elsewhere. It is also 
important to emphasize that dynamical inference of more 
sophisticated multi-dimensional models of the type [7(j 
can be addressed only in the frame of full Bayesian infer- 
ence of hidden dynamical variables. 

VII. CONCLUSION 

In the present paper we have introduced a technique for 
nonlinear dynamical inference of cardiovascular interac- 
tions from blood pressure time-series data. The method 
is applied to the simultaneous estimation of the dynami- 
cal couplings and noise strengths in a model of the nonlin- 
ear cardio-respiratory interaction. We have identified a 



simple nonlinear stochastic dynamical model of the car- 
diorespiratory interaction that describes, in framework 
of inverse modelling, the time-series data in a particu- 
lar frequency band. The method was validated using 
surrogate data obtained by numerically integrating the 
inferred model itself. We showed that main source of 
errors in the method is the decomposition of the blood 
pressure signal into two oscillatory components. We il- 
lustrate in the discussion that the dynamical model of 
the cardiorespiratory interaction identified in the present 
research can be related to the well-know beat-to-beat 
model of the cardiovascular control introduced by De- 
Boer and co-workers |59). The method introduced in 
this paper can be used to infer parameters of stochastic 
nonlinear dynamical models from observed phenomena 
across many scientific disciplines. 
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